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Past effort in computational techniques in internal flow systems has 
been concentrated on two-variable problems. This paper establishes a nu- 
merical method for the solution of three-variable problems and is applied 
here to rotational flows through ducts of various cross sections. 

An iterative scheme is developed, the main feature of which is the addi- 
tion of a duplicate variable to the forward component of velocity. Two 
forward components of velocity result from integrating two sets of first- 
order ordinary differential equations for the streamline curvatures, in 
intersecting directions across the duct. Two pseudo-continuity equations 
are introduced with source/sink terms, whose strengths are dependent on 
the difference between the forward components of velocity. When con- 
vergence is obtained, the two forward components of velocity are identi- 
cal, the source/sink terms are zero, and the original equations are 
satisfied. 

A computer program solves the exact equations and boundary condi- 
tions numerically. The method is economical and compares successfully 
with experiments on bent ducts of circular and rectangular cross section 
where secondary flows are caused by gradients of total pressure upstream. 


The presence of secondary-flow losses is well known. When a shear flow 
passes through a bend with a vorticity component directed toward the 
center of curvature, a secondary flow exists, transverse to the mean flow. 
The vorticity is produced by a velocity gradient in the flow approaching 
the bend. This velocity gradient may be produced by viscous losses up- 
stream and by nonuniform work being done on the fluid. The losses in a 
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secondary flow are due to the energy contained in the transverse flow, 
much of which is not recovered. (The presence of secondary flow may also 
cause subsequent parts entered by that fluid to run partly “off design.”) 

This attack on the secondary-flow problem solves the fully three- 
dimensional flow equations. The equations of three-dimensional fluid 
flow are intractable to analytic solution, even with the in viscid and steady 
flow assumptions. Until now, they have defied numerical solution due to 
the insufficient core size and speed of the past-generation computer and 
the lack of a numerical technique. The development of this three-dimen- 
sional method of solution was stimulated by the success of various 
two-dimensional numerical methods. The method is an extension of the 
two-variable streamline curvature method (refs. 1 and 2) . 

Although the method as presented is restricted to enclosed ducts, it is 
also possible to include repeat boundary conditions, thus enabling solu- 
tions of the turbomachinery blade passage flow to be obtained. 


DEVELOPMENT OF THE THREE-VARIABLE METHOD 

To economize on time and effort during the initial development of a 
three-variable method and to facilitate a clear understanding of the 
mechanisms involved, attention was restricted to incompressible flows 
and a simple geometry, for which experimental data was available (ref. 3) . 
The geometry is shown in figure 1. It consists of a rectangular duct which 
turns through any number of degrees on constant mean radius R m . 
Coordinates x and z are fixed in each plane of cross section and y is meas- 
ured along the centerline. 

An Eulerian approach to the equations is used, since the Lagrangian 
method, which is used in two-variable streamline curvature methods, is 
excessively complicated in three variables. It requires the storage and 
manipulation of expressions for two interacting families of stream surfaces 
and their interaction with the boundaries. 


Basic Equations 

Continuity 

1 d \( JL.\ ] 1 dv dw 

1 +x/R m dx \}+R m ) U yl-\-x/R m dy+dz 

Momentum 
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Figure 1. — Geometry of the duct. 
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These equations may be verified by considering a simple change of 
independent variables, r = R m +x and <t> = y/R m , which yields the well- 
known inviscid fluid flow equations in cylindrical polar coordinates. The 
components u, v, and w are physical velocity components normal to the 
local coordinate surfaces, and the static pressure p is understood to 
include the specific volume 1/p, which is constant. The term 1 /i? m is zero 
outside the bend region, corresponding to infinite radius of curvature. 


Manipulation of Equations 


Following usual streamline curvature procedure, the u and iv velocity 
components are replaced by new dependent variables X and p 


where 




u = hv 

(5) 

and 




W = y.V 

(6) 


(Most authors use tan X and tan p, but this is not necessary for the present 
analysis.) The “in-plane” components u and w, which are expected to be 
relatively small in an enclosed duct, are expressed as fractions of the 
dominating velocity component v, normal to the planes of cross section. 
No approximation is here implied, but the transformations (5) and (6) 
are singular when v is zero, a condition which must be avoided. Using 
(5) and (6), equations (1) through (4) become 
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Bernoulli’s equation is also derived from equations (8) through (10). 


dP 1 dP dP 
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(ID 


where the total pressure 

P = P+^ 2 (1 + X 2 +m 2 ) 


( 12 ) 


The five equations (8) through (12) are not independent, as Bernoulli’s 
equation is linearly dependent on the three momentum equations. One 
equation must be omitted, and (9) is selected since it is identical with 
(11) in the trivial case X = y = 0. 

Still adhering to the two-variable streamline curvature method, the 
dv/dy terms are eliminated from equations (8) and (10) with the aid of 
(7) to obtain 
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However, p is related to v 2 and P by (12). Using this equation in (13) 
and (14) 
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and 
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There are now four equations, (7), (11), (15) and (16), for the four 
dependent variables X, y, v, and P. Boundary conditions are required to 
close the system. 


Boundary Conditions 

It is necessary to appeal to the physics of the problem to obtain the 
correct boundary conditions. Some of these are obvious: X = 0 at the walls 
given by x= constant and y — 0 at the walls given by 2 = constant (the 
no-flow conditions). Upstream conditions are easily come by: straight 
shear flow where the static pressure is constant and either v or P is speci- 
fied at the inlet cross section. However, the conditions downstream are 
not so evident, being complicated by the presence of secondary flow. Two 
different downstream boundary conditions have been tried, both of which 
are sufficient to close the system of equations and boundary conditions 
from a numerical or computational point of view. The first condition is 
d\/dy = dy/dy = 0; the second is dv/dy = 0. The latter is a little more 
symmetric and converges faster, but both produce near-identical flow 
fields except over the last few computing planes. If, far downstream, there 
is a uniform swirling flow pattern, repeated at all subsequent planes of 
cross section, both boundary conditions are correct. 

Method of Solution of the Equations 

The extent to which the two-variable procedure may be followed has 
now been reached. Examination of the equations indicates the following: 

( 1 ) Bernoulli’s equation (eq. (11)). Given values for X and y through- 
out the flow field, P may be calculated from the starting values at the 
inlet cross section. 

(2) The momentum equations (eqs. (15) and (16)). Either of these 
may be integrated for v when X, y, and their derivatives are known. 

(3) The continuity equation (eq. (7)). Assuming that v is given 
throughout the flow field, this equation may be integrated for X if y is 
known or for y if X is known. 

These integrations will be for linear, first-order, ordinary differential 
equations with nonconstant coefficients. Bernoulli’s equation is written 
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DP/Dt = 0 along a streamline to fit this classification. It is inconvenient, 
though, to have two equations, either of which may be solved explicitly 
for v, and only one equation to solve for both X and /x. Two methods of 
solution have been tried. 

First Method of Solution 

A pattern similar to that proposed by Wu (ref. 4) was followed. 
Separate sets of two-dimensional solutions were sought , with an iterative 
procedure connecting them as shown in figure 2. This did not work. 
Alternating between one set of two-dimensional solutions and the other 
is not sufficient to produce convergence since neither solution “realizes” 
that it is not the same as the other. The information conveyed between 
the solutions is not sufficient to produce convergence. 

Although Wu’s proposals differ in that his two-dimensional solutions 
are calculated on Si and St stream surfaces using two stream functions, 
the method follows the pattern suggested by Wu. A few variations on this 
method have also been tried, but without success. This suggests that for a 
method to have any chance of success it must “know” about the “error” 
or difference between separate two-dimensional solutions, v x — v z for 
example, and act on this information until the error is reduced to zero. 

Second Method of Solution 

Let the result of integrating equation (15) in the ^-direction be v x , 
and the result of integrating equation (16) in the z-direction be v z . The 
error v x — v‘ is related to a static pressure difference by equation (12). 
Physically, this pressure difference will change the curvature of the 
streamlines, and thus X and fi must be influenced by v x —v z . The best choice 
seems to be the replacement of the continuity equation (eq. (7)) by the 
two equations, 

d 

fix 

and 

fi 

fix 

where /? is a constant. In these equations, the right-hand sides represent 
source/sink terms and each reduces to the continuity equation (eq. (7)) 
when v x = v x . One additional equation and one additional unknown have 
been introduced and now (17) is integrated directly for X, (18) for n, 
(15) for v x , (16) for v z , and (11) for P. In (17) and (18), v x and v z are 
selected appropriately to make the boundary conditions for the velocity 
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integrations explicit. The essential feature of these equations is that when 
v x = v z they reduce to the physically correct equations. 


START 



END 


Figure 2. — Iterative procedure for Method 1. 
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Revised Boundary Conditions 

Although the boundary conditions described are sufficient to determine 
a solution, they do not lend themselves to an easy integration of the 
equations in their present form. First, (17) and (18) are first-order 
equations; each requires only one boundary condition for X and y and 
there are two for each. Second, (15) and (16) also require one boundary 
condition each for v x and v z and there are none. This matter is rectified by 
the requirements that two boundary conditions may be satisfied by each 
of the first-order pseudo-continuity equations (eqs. (17) and (18)). 
These requirements are found as follows. Equations (17) and (18) are 
integrated, first with respect to x from x = — \X across the duct to x = 5 X 
and then with respect to y from y = 0 at the inlet to some station y , 
to yield 


A/ZX 

I v 1 dx-\- 

•'-1/2X 



(yv z ) dxdy = 0 


(19) 


and 


A/iX r v A/ix / r \ Q 

[ v*dx+ [ f (l+— )- (/»*) dxdy = 0 (20) 

—1/2X *'0 •'-1/2X V d z 


where the p(v z —v*) terms have been omitted, and the boundary condi- 
tions X = 0 at x= ±%X have been incorporated. Alternatively, a repeat 
condition, \(-\X,y,z) = X(§X,J/,z) and v*(-\X,y,z) = v x (\X,y,z), yields 
the same results. The requirement (19) is used as a boundary condition 
for equation (15). The procedure is repeated with the roles of x and 2 
interchanged to obtain a similar requirement for the other pseudo- 
continuity equation and a boundary condition for (16). 

In general, sets of coupled partial differential equations cannot be put 
into explicit form, so it is necessary to select one variable in an equation 
and guess or assume values for all others. Each of the variables must 
take its turn as the unknown in one of the equations. When all variables 
have been found, the equations are solved again and this iterative pro- 
cedure is continued until convergence is obtained. 


Iterative Scheme 

Figure 3 shows the iterative scheme. Each block represents the integra- 
tion of the appropriate equation for the unknown variable throughout 
the entire flow field. In each integration, the most up-to-date values are 
used for all other variables. This scheme is chosen for its simplicity and 
because it also simplifies the boundary conditions of (19) and (20). 
Between the calculation of y and the next calculation of X, v* remains 
unchanged; hence, on subtracting equations (19) and (20) 
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AI1X .1/2X 

/ v x dx = v‘ dx (21) 
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This is the boundary condition used for v x . Similarly, the boundary 
condition for v‘ is 



The right-hand sides of these 
culations. 
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Figure 3. — Iterative procedure for 
Method 8. 
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Uniqueness 

These boundary conditions for the velocity are not of a familiar type, 
and the existence of a unique solution must be explored. Equation (15) 
is like the equation 

— v i = v i g{x) +/(x) (23) 

dx 

A solution is sought, subject to the condition 

f vdx = Q (24) 

•'o 


If two solutions, v\ and v^, exist then 

4- (vS-vS) = (iY 2 -f 2 2 ) g(x) 
dx 


subject to 



The solution of (25) is 


vi 2 — 1 ’ 2 2 = E exp 



(25) 


(26) 


(27) 


where E is determined by (26) as follows 


E 


L 


exp (j g{s) dsj 


V\ + l'2 


dx = 0 


(28) 


The further restriction that v>0 is necessary for uniqueness. Now, the 
integrand of (28) is always positive; whence E = 0 and *<’i = in- 


stability 

Theoretically, it is only possible to perform a stability analysis for 
trivial flows where the total pressure is constant, but the resulting criterion 
is found to have general application. A straight flow without shear is 
considered, where v x =v z =V and X = n = 0 is the required solution. Small 
perturbations from this trivial solution are examined and the following 
stability criterion is obtained 



1+2 ( ara [! (° 63 0 ^] 


<1 


(29) 
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where X is the duct width, Ay is the grid spacing in the y direction, and 
r is the relaxation factor on the velocity. 

A computer program has been written in FORTRAN IV to solve the 
problem as outlined. An experiment with a 5-inch by 5-inch 90° bend on a 
15-inch mean radius, where the velocity varied across the duct from about 
30 feet/seconu to 80 feet/seeond at inlet, was simulated. The results, in 
the form of P 0 contours, are presented for comparison with experiment. 
The small discrepancies can be accounted for by the presence and separa- 
tion of the boundary layer. Otherwise, an adequate prediction of the flow 
is obtained. 


Comparison With Experiment 

Numerical computations have been carried out for the experiments of 
Joy (ref. 3) for ducts of rectangular cross section bending through 90° 
and subject to substantial inlet total head variations across the duct. 
A comparison has also been made with ducts of circular cross section 
(Eichenburger reported in ref. 5). The theory presented in this paper is 
directly applicable to the rectangular duct but requires modification to 
the circular geometry although the equations are of a similar form. 

In figure 4 the total pressure contours at inlet to the duct are presented. 
In figures 5, 6, and 7 the computed contours are compared with experi- 
ment at three stations down the duct. 

Station 1 6 inches upstream of the bend 

Station 2 30° of turning 

Station 3 60° of turning 

Station 4 90° of turning 

The duct is 5 inches by 10 inches in cross section with a mean radius of 
15 inches. For consistency with Joy, the total pressure contours are 
labelled as velocity contours computed on the assumption of constant 
static pressure. Similar comparisons are shown for the circular-cross- 
sectioned duct in figures 8, 9, and 10. The duct is of 6 inches diameter and 
30 inches mean radius. 

In general, the experimental contours are predicted by the theory. 
For the circular duct, the agreement is particularly good except in the 
immediate vicinity of the wall where the viscous forces in the boundary 
layer are dominant, causing reductions in total pressure. The discrepancies 
in predictions for the rectangular duct near the inside of the bend are 
probably due to the occurrence of separation of the boundary layer near 
Station 3. 

A measure of the convergence of the numerical procedure for the 
rectangular bend is presented in figure 11 showing good convergence 
after 58 cycles. This procedure took 14 minutes on an IBM 360/65. 
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LIST OF SYMBOLS 

p Static pressure 

P 0 Total or stagnation pressure 
u, v, w Velocity components 
x, y, z Coordinates 

d Strength of source/sink distribution 

X, pi Flow directions as defined 
p Density 


15* x 90° BEND 

STATION I - LOOKING UPSTREAM 



Figure 4. — Velocity contours in rec- 
tangular duct. 
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15* X 90° BEND 

STATION 2 - LOOKING UPSTREAM 



EXPERIMENT THEORY 

Figure 5. — Velocity contours in rectangular duct; comparison between theory and 

experiment. 


15* x 90° BEND 

STA TI ON 3 -LOOKING UPSTREAM 



EXPERIMENT THEORY 

Figure 6. — Velocity contours in rectangular duct; comparison between theory and 

experiment. 
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15* x 90° BEND 

STATION 4 - LOOKING UPSTREAM 



Figure 7. — Velocity contours in rectangular dud; comparison between theory and 

experiment. 
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Total pressure Contours, 
due to EICHENBURGER, in a 
6-inch-diameter circular duct, 
bent on a mean radius of 
30 inches. 


Figure 8. — Total pressure contours in circular duct. 
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30 ° 



Figure 9. — Total pressure contours in circular duct; comparison between theory and 

experiment. 



Figure 10. — Total pressure contours in circular duct; comparison between theory and 

experiment. 
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Figure 11. — Convergence of 
computing scheme. 
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DISCUSSION 


H. MARSH (Cambridge University) : The authors have successfully 
solved a major problem in the calculation of internal flows, namely the 
three-dimensional flow in a duct. Until recently, it has not been possible 
to solve this problem due to the lack of adequate computing facilities. It 
would be interesting to know the computer storage requirements for this 
program. 

Two methods of solution have been examined, but only the successful 
second method is described in detail. The first method is similar to that 
proposed by Wu (ref. 4). Until now, this has been considered a viable 
method for calculating the full three-dimensional flow field. The authors 
have investigated this technique and have found that they could not 
obtain convergence. This negative result is extremely important and it 
deserves a more detailed discussion. Smith (ref. D-l) has described the 
methods which are available for calculating the two separate two-dimen- 
sional flow fields. Until now, it has been assumed that by alternating 
between the two solutions, the full three-dimensional flow field might be 
calculated. It would be helpful if the authors would give more details of 
the basis for their conclusion that the first method of solution does not 
work. 

In the second method of solution, the error v z — v‘ is related to a static 
pressure difference, but it is not clear why this term should have any 
physical significance. The replacement of the continuity equation by two 
equations with source/sink terms is a numerical technique which is used 
in order to obtain a convergent solution. It is therefore unlikely that the 
intermediate values of the error v z —v z have any physical meaning. 

In the derivation of the boundary conditions, the authors have omitted 
the source/sink term but have not discussed this point. Perhaps they 
w’ould outline their argument for neglecting these terms. It is possible to 
argue that any convenient boundary condition can be used, provided that 
it approaches the true boundary condition as the solution converges. 

This is a major contribution to methods of flow calculation and the 
authors must be congratulated on their presentation in this paper. If this 
work can be extended to include compressibility, then it would provide a 
single comprehensive technique for calculating inviscid three-dimen- 
sional duct flows. 

W. R. HAWTHORNE (Cambridge University) : I agree with what 
Mr. Stuart says. I think the work of Rowe (ref. D-2) should be referred 
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to in this excellent paper which seems to me to be a substantial break- 
through on three-dimensional and secondary flow calculations. But I also 
want to raise the question of getting the right downstream boundary 
conditions. It isn’t clear from what the author was saying how and where 
the boundary conditions were established. Were they established at 30°, 
40°, or down the bend? In one case, he showed a section 1 foot downstream 
from the bend. How far downstream could you go before you get the 
right downstream conditions? 

STUART and HETHERINGTON (authors) : In reply to H. Marsh, 
we would like to state the following : 

(1) The computations were performed on an IBM 360/65 and re- 
quired between 120 000 and 180 000 bytes (depending on whether the 
program was overlayed or not). Typical execution times were 10 to 14 
minutes CPU. 

(2) The conclusion that simple alternation between two separate 
two-dimensional solutions does not produce the required three-dimen- 
sional flow field is based on our failure to make such methods produce 
identical fields for the axial velocity from both two-dimensional solutions, 
in the absence of the /3(v x — v z ) terms in equations (17) and (18). For test 
computations with (3 = 0, convergence has not been obtained, and over a 
considerable portion of the flow field (about half), near the start of the 
bend, the secondary flow turned in the wrong direction. This even 
propagated upstream where no secondary flow is to be expected. 

(3) We agree that physically no terms 0(v z —v z ) exist. The argument 
for the use of such a term as a numerical device is as follows : Physically, 
fluid will tend to flow from high pressure toward low pressure regions until 
the pressure gradient is balanced by acceleration. Now, an imbalance may 
exist between the pressure gradient of one two-dimensional solution and 
the acceleration or curvature terms of the other (since these balance their 
own pressure gradient, which is not necessarily identical to that of the 
first solution). This imbalance between respective pressure gradients is 
related to v x —v z , which term or “error” is used to change X and y ac- 
cordingly in the psendo-continuity equations (17) and (18). 

(4) The neglect of some terms 0(v x —v z ) in the derivation of the 
revised boundary conditions is justified as follows: In practice, the pseudo- 
continuity equations (17) and (18) arc solved with additional source and 
sink terms S x (y,z) and S z (x,y ) to allow for the effect of the terms omitted. 
Now, it must be shown that these terms vanish when convergence is 
obtained. If the two solutions v x and v z converge (remember that these are 
both the axial velocity, and not components in the x and z directions), 
then subtracting equation (17) from equation (18) yields, at most 


S x (y,z) =S z (x,y ) =f(y) 
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Integrating equation (17) or (18) over a plane of cross section yields 

f(y) j dxdz = ^ 

where Q is the mass flow. Since Q is constant, f(y ) = 0 = S x = S z . 

(5) In its present form, the method has been extended to compressible 
flows (Mach number less than 0.98) and the computer program can 
handle the following problems : 

(a) Bent rectangular ducts 

(b) Bent circular ducts, including S-bends 

(c) Three-dimensional flow through a turbomachinery blade 
passage specified by random points, including rotors. 

Arbitrary values for total temperature, total pressure, and static pressure, 
varying across the inlet section may be specified as input data to the 
program. Future work in the Department of Mathematics at the Uni- 
versity of Aston will attempt to extend the method further, to include 
viscous and turbulent flows. 

As regards the right boundary conditions mentioned by Sir William 
Hawthorne, we would like to offer the following reply. 

The computing mesh is usually extended two or three planes further 
downstream of the region of interest (there being expense involved in 
using too many) , but from the calculations we have done, the condition 
downstream does not have much effect two or three planes upstream 
(i.e., about three pipe diameters) of where the downstream boundary 
condition is applied. The flow fields for the two different downstream 
boundary conditions described were within 0. 1 percent of being identical 
two planes upstream of where the conditions were established. 
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